from numpy.lib.shape_base import _replace_zero_by_x_arrays
from sympy import *

x=Symbol('x')
y=Symbol('y')

print("derivative of linear basis function:")
p1=-x-y+1
p2=x
p3=y
print("\np1:")
print(simplify(p1.diff(x,1,y,0)))
print(simplify(p1.diff(x,0,y,1)))
print(simplify(p1.diff(x,1,y,1)))

print("\np2:")
print(simplify(p2.diff(x,1,y,0)))
print(simplify(p2.diff(x,0,y,1)))
print(simplify(p2.diff(x,1,y,1)))

print("\np3:")
print(simplify(p3.diff(x,1,y,0)))
print(simplify(p3.diff(x,0,y,1)))
print(simplify(p3.diff(x,1,y,1)))

x1=Symbol('x1')
x2=Symbol('x2')
x3=Symbol('x3')
y1=Symbol('y1')
y2=Symbol('y2')
y3=Symbol('y3')
x4=Symbol('x4')
y4=Symbol('y4')

xh=Symbol('xh')
yh=Symbol('yh')


f1=x-((x2-x1)*xh+(x3-x1)*yh+x1)
f2=y-((y2-y1)*xh+(y3-y1)*yh+y1)
r=solve([f1,f2],[xh,yh])
print(simplify(r[xh]))
print(simplify(r[yh]))

print("2D quadratic basis function")

p1=2*x**2+2*y**2+4*x*y-3*y-3*x+1
p2=2*x**2-x
p3=2*y**2-y
p4=-4*x**2-4*x*y+4*x
p5=4*x*y
p6=-4*y**2-4*x*y+4*y

print("p1:")

print("p1/x={}".format(p1.diff(x,1,y,0)))
print("p1/y={}".format(p1.diff(x,0,y,1)))

print("p1/x.p1/y={}".format(p1.diff(x,1,y,1)))
print("p1^2/x^2={}".format(p1.diff(x,2,y,0)))
print("p1^2/y^2={}".format(p1.diff(x,0,y,2)))

print("p1^3/x^3={}".format(p1.diff(x,3,y,0)))
print("p1^3/y^3={}".format(p1.diff(x,0,y,3)))
print("p1/x.p1^2/y^2={}".format(p1.diff(x,1,y,2)))
print("p1^2/x^2.p1/y={}".format(p1.diff(x,2,y,1)))

print("p2:")

print("p2/x={}".format(p2.diff(x,1,y,0)))
print("p2/y={}".format(p2.diff(x,0,y,1)))

print("p2/x.p2/y={}".format(p2.diff(x,1,y,1)))
print("p2^2/x^2={}".format(p2.diff(x,2,y,0)))
print("p2^2/y^2={}".format(p2.diff(x,0,y,2)))

print("p2^3/x^3={}".format(p2.diff(x,3,y,0)))
print("p2^3/y^3={}".format(p2.diff(x,0,y,3)))
print("p2/x.p2^2/y^2={}".format(p2.diff(x,1,y,2)))
print("p2^2/x^2.p2/y={}".format(p2.diff(x,2,y,1)))

print("p3:")

print("p3/x={}".format(p3.diff(x,1,y,0)))
print("p3/y={}".format(p3.diff(x,0,y,1)))

print("p3/x.p3/y={}".format(p3.diff(x,1,y,1)))
print("p3^2/x^2={}".format(p3.diff(x,2,y,0)))
print("p3^2/y^2={}".format(p3.diff(x,0,y,2)))

print("p3^3/x^3={}".format(p3.diff(x,3,y,0)))
print("p3^3/y^3={}".format(p3.diff(x,0,y,3)))
print("p3/x.p3^2/y^2={}".format(p3.diff(x,1,y,2)))
print("p3^2/x^2.p3/y={}".format(p3.diff(x,2,y,1)))


print("p4:")

print("p4/x={}".format(p4.diff(x,1,y,0)))
print("p4/y={}".format(p4.diff(x,0,y,1)))

print("p4/x.p4/y={}".format(p4.diff(x,1,y,1)))
print("p4^2/x^2={}".format(p4.diff(x,2,y,0)))
print("p4^2/y^2={}".format(p4.diff(x,0,y,2)))

print("p4^3/x^3={}".format(p4.diff(x,3,y,0)))
print("p4^3/y^3={}".format(p4.diff(x,0,y,3)))
print("p4/x.p4^2/y^2={}".format(p4.diff(x,1,y,2)))
print("p4^2/x^2.p4/y={}".format(p4.diff(x,2,y,1)))

print("p5:")

print("p5/x={}".format(p5.diff(x,1,y,0)))
print("p5/y={}".format(p5.diff(x,0,y,1)))

print("p5/x.p5/y={}".format(p5.diff(x,1,y,1)))
print("p5^2/x^2={}".format(p5.diff(x,2,y,0)))
print("p5^2/y^2={}".format(p5.diff(x,0,y,2)))

print("p5^3/x^3={}".format(p5.diff(x,3,y,0)))
print("p5^3/y^3={}".format(p5.diff(x,0,y,3)))
print("p5/x.p5^2/y^2={}".format(p5.diff(x,1,y,2)))
print("p5^2/x^2.p5/y={}".format(p5.diff(x,2,y,1)))

print("p6:")

print("p6/x={}".format(p6.diff(x,1,y,0)))
print("p6/y={}".format(p6.diff(x,0,y,1)))

print("p6/x.p6/y={}".format(p6.diff(x,1,y,1)))
print("p6^2/x^2={}".format(p6.diff(x,2,y,0)))
print("p6^2/y^2={}".format(p6.diff(x,0,y,2)))

print("p6^3/x^3={}".format(p6.diff(x,3,y,0)))
print("p6^3/y^3={}".format(p6.diff(x,0,y,3)))
print("p6/x.p6^2/y^2={}".format(p6.diff(x,1,y,2)))
print("p6^2/x^2.p6/y={}".format(p6.diff(x,2,y,1)))

print("2D bilinear basis functions:")

p1=(1-x-y+x*y)/4
p2=(1+x-y-x*y)/4
p3=(1+x+y+x*y)/4
p4=(1-x+y-x*y)/4

print("p1:")

print("p1/x={}".format(p1.diff(x,1,y,0)))
print("p1/y={}".format(p1.diff(x,0,y,1)))

print("p1/x.p1/y={}".format(p1.diff(x,1,y,1)))
print("p1^2/x^2={}".format(p1.diff(x,2,y,0)))
print("p1^2/y^2={}".format(p1.diff(x,0,y,2)))

print("p2:")

print("p2/x={}".format(p2.diff(x,1,y,0)))
print("p2/y={}".format(p2.diff(x,0,y,1)))

print("p2/x.p2/y={}".format(p2.diff(x,1,y,1)))
print("p2^2/x^2={}".format(p2.diff(x,2,y,0)))
print("p2^2/y^2={}".format(p2.diff(x,0,y,2)))


print("p3:")

print("p3/x={}".format(p3.diff(x,1,y,0)))
print("p3/y={}".format(p3.diff(x,0,y,1)))

print("p3/x.p3/y={}".format(p3.diff(x,1,y,1)))
print("p3^2/x^2={}".format(p3.diff(x,2,y,0)))
print("p3^2/y^2={}".format(p3.diff(x,0,y,2)))

print("p4:")

print("p4/x={}".format(p4.diff(x,1,y,0)))
print("p4/y={}".format(p4.diff(x,0,y,1)))

print("p4/x.p4/y={}".format(p4.diff(x,1,y,1)))
print("p4^2/x^2={}".format(p4.diff(x,2,y,0)))
print("p4^2/y^2={}".format(p4.diff(x,0,y,2)))

print("affine mapping for 2D bilinear element:")

a=Symbol('a')
b=Symbol('b')
c=Symbol('c')
d=Symbol('d')
e=Symbol('e')
f=Symbol('f')

x1=1
y1=-1
x2=1
y2=-0.9
x3=0.9
y3=-0.9
x4=0.9
y4=-1.0

f1=a*x1+b*y1+c+1
f2=a*x2+b*y2+c-1
f3=a*x3+b*y3+c-1

r=solve([f1,f2,f3],[a,b,c])
print("xh:")
print(r)
print("x4:")
print(r[a]*x4+r[b]*y4+r[c])

f4=d*x1+e*y1+f+1
f5=d*x2+e*y2+f+1
f6=d*x3+e*y3+f-1

r=solve([f4,f5,f6],[d,e,f])
print("yh:")
print(r)

print("y4")
print(r[d]*x4+r[e]*y4+r[f])

print("x,y")
f7=xh-(a*x+b*y+c)
f8=yh-(d*x+e*y+f)

r=solve([f7,f8],[x,y])
print(r)



print("dx/dxh")
print(r[x].diff(xh,1,yh,0))
print("dx/dyh")
print(r[x].diff(xh,0,yh,1))
print("dy/dxh")
print(r[y].diff(xh,1,yh,0))
print("dy/dyh")
print(r[y].diff(xh,0,yh,1))
print("dydy")
print(r[y].diff(xh,1,yh,1))


f1=x-(-b*yh+e*xh)/(a*e-b*d)
f2=y-(a*yh-d*xh)/(a*e-b*d)
r=solve([f1,f2],[xh,yh])
print(r)

print("2D bilinear")
h1=Symbol('h1')
h2=Symbol('h2')
f1=xh-(2*x-2*x1-h1)/h1
f2=yh-(2*y-2*y1-h2)/h2
r=solve([f1,f2],[x,y])
print(r)















